clear all

load ParEst_eh0_time0_Mex_guys ParEst_eh0_time0_Mex_guys
load ParEst_eh0_time2_Mex_guys ParEst_eh0_time2_Mex_guys
load MWP_a
load stock_time2_age0_eh0_comp1.out
load W_time2_age0_eh0_comp11.out
load W_time2_age0_eh0_comp1.out

Par1 = zeros(14,1);
Par1(1:13) = ParEst_eh0_time0_Mex_guys;
Par1(11) = 10^Par1(11);
Par1(12:13)=1000*Par1(12:13);
Par2 = ParEst_eh0_time2_Mex_guys;
Par1(14)=500*Par2(1);
Par1(15:18)=Par2(2:5);

% Parameters labels - men
df_1=Par1(1);
di_1=Par1(2);
lnf_1=Par1(3);
lni_1=Par1(4);
lfi_1=Par1(5);
lif_1=Par1(6);
alphaf=Par1(7);
betaf=Par1(8);
alphai=Par1(9);
betai=Par1(10);
theta=Par1(11);
a=Par1(12);
b1 = Par1(13);
gamma=Par1(14); 
alphaf_2=Par1(15);
betaf_2=Par1(16);
alphai_2=Par1(17);
betai_2=Par1(18);

w_f_g=W_time2_age0_eh0_comp11; % Read wages, g and f distributions
wf_1=w_f_g(:,1);        % Wage in formal sector
wi_1=w_f_g(:,2);        % Wage in informal sector
wf_2=w_f_g(:,7);      % Wage in formal sector
wi_2=w_f_g(:,8);      % Wage in informal sector

w_mom=W_time2_age0_eh0_comp1; % Read wage moments
meanwi_1=w_mom(1); meanwf_1=w_mom(2);   meanwi_2=w_mom(3);  meanwf_2=w_mom(4);  sdwi_1=w_mom(5); sdwf_1=w_mom(6);  sdwi_2=w_mom(7);   sdwf_2=w_mom(8);
p10wi_1=w_mom(9);  p10wf_1=w_mom(10); p10wi_2=w_mom(11); p10wf_2=w_mom(12); p25wi_1=w_mom(13); p25wf_1=w_mom(14); p25wi_2=w_mom(15); p25wf_2=w_mom(16); 
p50wi_1=w_mom(17); p50wf_1=w_mom(18); p50wi_2=w_mom(19); p50wf_2=w_mom(20); p75wi_1=w_mom(21); p75wf_1=w_mom(22); p75wi_2=w_mom(23); p75wf_2=w_mom(24); 
p90wi_1=w_mom(25); p90wf_1=w_mom(26); p90wi_2=w_mom(27); p90wf_2=w_mom(28); minwi_1=w_mom(29); minwf_1=w_mom(30); minwi_2=w_mom(31); minwf_2=w_mom(32); 
maxwi_1=w_mom(33); maxwf_1=w_mom(34); maxwi_2=w_mom(35); maxwf_2=w_mom(36);
wage_moments_data=[meanwi_1; meanwf_1; sdwi_1; sdwf_1]; % guys

stock_mom=stock_time2_age0_eh0_comp1; % Read stocks
stock_ANT=stock_mom'; 
stock = [stock_ANT(1) + stock_ANT(2) + stock_ANT(3);...
    stock_ANT(4) + stock_ANT(6) + stock_ANT(7); ...
    stock_ANT(5) + stock_ANT(8) + stock_ANT(9)]; % guys

HHstates_data=stock;

n=50;
 % solve dist of wage offers, value functions and simulate trajectory
[Ff_1,ff_1,Fi_1,fi_1]=solve_wage_offer_dist(Par1,wf_1,wi_1,n);
[Wf,Wi,Wn,iterVF]=solve_value_function_A1(Par1,Ff_1,wf_1,Fi_1,wi_1,n);
[HHstates,wage_moments,transitions]=data_simulation_A12_end(Par1,ff_1,wf_1,fi_1,...
    wi_1,Wf,Wi,Wn);


HH_avg_Inc = p50wi_1; % median wage of informal men
MWP_gamma= gamma/exp(-theta*HH_avg_Inc);
save MWP_gamma MWP_gamma

%Preference parameters
Preference = {'theta';'b1';'a';'MWP a';'gamma';'MWP gamma'};
Estimates = [theta; b1; a; MWP_a;gamma; MWP_gamma];
Preferences_table = table(Preference, Estimates);
writetable(Preferences_table, 'Preferences_table.xls')

%Wage offer parameters 1st stage
Wage_offer = {'alpha_formal'; 'beta_formal';'alpha_informal';'beta_informal'};
Estimates = [alphaf; betaf; alphai; betai];
Wage_offer_par_table = table(Wage_offer, Estimates);
writetable(Wage_offer_par_table, 'Wage_offer_par_table.xls')

% Generate mean LOG WAGES - F distribution 
mean_wf_1f=log(sum(wf_1.*ff_1));
mean_wi_1f=log(sum(wi_1.*fi_1));
mean_wage_offer_time1=[mean_wf_1f;mean_wi_1f];
save mean_wage_offer_time1 mean_wage_offer_time1

%Wage offer parameters 2nd stage
Wage_offer = {'alpha_formal'; 'beta_formal';'alpha_informal';'beta_informal';'';'mean_wf_1f';'mean_wi_1f'};
Estimates = [alphaf_2; betaf_2; alphai_2; betai_2;NaN;mean_wage_offer_time1];
Wage_offer_par_table = table(Wage_offer, Estimates);
writetable(Wage_offer_par_table, 'Wage_offer_par_table_2nd_stage.xls')

x=wage_moments_data;
wage_moments_data1=[x(2);x(4);NaN;NaN;...	
x(1);x(3)];
w=wage_moments(1:4);
wage_moments1=[w(2);w(4);NaN;NaN;... 
w(1);w(3)];
% Model fit (for paper)
Model_fit_labels = {'m_f';'m_i';'m_n';'';'';'meanwf_1';'sdwf_1';'';'';'meanwi_1';'sdwi_1'};
ModelFit_actual= [HHstates_data;NaN;NaN;log(wage_moments_data1)];
ModelFit_model= [HHstates;NaN;NaN;log(wage_moments1)];
Model_fit_table = table(Model_fit_labels,ModelFit_actual,ModelFit_model);
writetable(Model_fit_table, 'Model_fit_table_paper.xls')